#!/bin/bash
# ----------------------------------------------------------------
#   3-B0_PrepareFieldmap.sh: Prepare B0 field map file, which will
#	be used for un-warping later on
#
#	Usage:
#		3-B0_PrepareFieldmap.sh B0_phase1 B0_phase2 B0_mag
#   initializes the Adaptation fMRI experiment. This experiment 
#   effect of stimulus history on the current evoked response,
#   on serial position -1 to -N are manipulated in a 
#
#
#   Inputs
#   -------
#   B0_phase1, B0_phase2, B0_mag: phase and magnitude images
#
#
#   Examples
#   -------
#   3-B0_PrepareFieldmap.sh B0_phase1 B0_phase2 B0_mag
#
# ----------------------------------------------------------------
# ``Redistribution and use, with or without modification, are
#	permitted provided that the list of contributors below is
#	preserved.``
#
#   Original code by Marcelo G Mattar (09/04/2012)
#	mattar@sas.upenn.edu
# ----------------------------------------------------------------


# --- Check input command ---
if [ $# -ne 3 ]; then
cat << EOF
USAGE: `basename $0` B0_phase1 B0_phase2 B0_mag
EOF
	exit 1
fi

# --- Read input parameters ---
B0_phase1File=$1
B0_phase2File=$2
B0_magFile=$3


# --- Define some basic parameters ---
DELTA_TE=2.58
SCALING=4096
B0DIR=`dirname $B0_magFile`
FIELDMAPDIR=$B0DIR/fieldmap
BET_THRESHOLD=0.7

# --- Keep a copy of calling directory, so that relative pathnames are properly used
OCD=$PWD


# -----------------------------------------------------
# 1) --- CREATE FIELDMAP DIRECTORY ---
# -----------------------------------------------------

# Create /FieldMap directory
if [ -d "$FIELDMAPDIR" ]; then
	echo
	echo "WARNING: FIELDMAPDIR directory $FIELDMAPDIR already exists!!!"
	read -p "Overwrite (y/n)? " -n 1
	if [[ ! $REPLY =~ ^[Yy]$ ]]
	then
		echo "Exiting..."
		echo
		exit 1
	fi
else
	echo
	mkdir $FIELDMAPDIR
fi

echo


# -----------------------------------------------------
# 2) --- EXTRACT THE BRAIN OF THE MAGNITUDE IMAGE ---
# -----------------------------------------------------

echo "Performing brain extraction on magnitude image..."

bet $B0_magFile $FIELDMAPDIR/B0_mag_BET_0.5.nii.gz -R -f 0.5 -g 0 -o
bet $B0_magFile $FIELDMAPDIR/B0_mag_BET_0.6.nii.gz -R -f 0.6 -g 0 -o
bet $B0_magFile $FIELDMAPDIR/B0_mag_BET_0.7.nii.gz -R -f 0.7 -g 0 -o

echo "Done!"
echo


# -----------------------------------------------------
# 3) --- CONVERT PHASE IMAGES TO RADIANS ---
# -----------------------------------------------------

echo "Converting phase images to radians..."

# Apparently, the final range of values here must be between -pi and pi, instead of 0 and 2pi.

#fslmaths $B0_phase1File -mul 3.14159 -div $SCALING -add 3.14159 $FIELDMAPDIR/B0_phase1_rad_wrapped.nii.gz -odt float
#fslmaths $B0_phase2File -mul 3.14159 -div $SCALING -add 3.14159 $FIELDMAPDIR/B0_phase2_rad_wrapped.nii.gz -odt float

fslmaths $B0_phase1File -mul 3.14159 -div $SCALING $FIELDMAPDIR/B0_phase1_rad_wrapped.nii.gz -odt float
fslmaths $B0_phase2File -mul 3.14159 -div $SCALING $FIELDMAPDIR/B0_phase2_rad_wrapped.nii.gz -odt float

echo "Done!"
echo


# -----------------------------------------------------
# 4) --- UNWRAP PHASE IMAGES ---
# -----------------------------------------------------

echo "Unwrapping phase images..."

prelude -a $FIELDMAPDIR/B0_mag_BET_${BET_THRESHOLD}.nii.gz -p $FIELDMAPDIR/B0_phase1_rad_wrapped.nii.gz -o $FIELDMAPDIR/B0_phase1_rad_unwrapped.nii.gz
prelude -a $FIELDMAPDIR/B0_mag_BET_${BET_THRESHOLD}.nii.gz -p $FIELDMAPDIR/B0_phase2_rad_wrapped.nii.gz -o $FIELDMAPDIR/B0_phase2_rad_unwrapped.nii.gz

echo "Done!"
echo


# -----------------------------------------------------
# 5) --- CALCULATE FIELDMAPS IN RADIANS ---
# -----------------------------------------------------

echo "Generating fieldmaps..."

# Apparently, the right thing to do is phase2 - phase1 (if phase encoding direction is -y)
# That's what seems to give us images that are similar to what FSL had in their tutorial

fslmaths $FIELDMAPDIR/B0_phase2_rad_unwrapped.nii.gz -sub $FIELDMAPDIR/B0_phase1_rad_unwrapped.nii.gz -mul 1000 -div $DELTA_TE $FIELDMAPDIR/B0fieldmap1_rads.nii.gz -odt float
fslmaths $FIELDMAPDIR/B0_phase1_rad_unwrapped.nii.gz -sub $FIELDMAPDIR/B0_phase2_rad_unwrapped.nii.gz -mul 1000 -div $DELTA_TE $FIELDMAPDIR/B0fieldmap2_rads.nii.gz -odt float

echo "Done!"
echo

echo "Fieldmaps saved at ${FIELDMAPDIR}"





cd $OCD

exit 0
